#' Total, Direct, and Indirect Spatiotemporal Effects for the STADL (sy0,ty1) Model.
#'
#' \code{st_effects} calculates both the short and long-run total, direct and indirect effects of shocks to sample units' covariates using the last cross-section of the data.
#'
#' @param spobj An object created by the \code{ntspreg} function for a STADL (sy0,ty1) Model.
#' @param wm a list created by the \code{make_ntspmat} function that includes the spatial weights    #' matrix.
#' @param tlag Name of a variable that identifies the temporal lag of an outcome.
#' @param covariate Name of a variable to be shocked.
#' @return The output will be a dataframe.
#'
#'
#' @import dplyr
#' @importFrom stats na.omit
#' @importFrom sf as_Spatial
#' @importFrom Matrix bdiag
#' @importFrom lubridate year
#' @import cshapes
#' @import stargazer
#' @export


#' @examples
#'  \dontrun{
#' fx <- st_effects(sar,wm,tlag,x)
#' }

st_effects_new <- function(spobj,wm,tlag,covariate) {

  call <- match.call()
  df <- wm[[1]]
  yi<- df[,wm[[4]]]
  yit <- as.character(wm[[4]])
  dim0 <- as.numeric(length(which(eval(parse(text = noquote(paste("df$",yit,sep="")))) < max(yi))))
  dim1 <- as.numeric(length(which(df$year == max(yi))))

  W <- as.matrix(wm[[2]]/rowSums(wm[[2]]))
  W <- W[(dim0+1):(dim0+dim1),(dim0+1):(dim0+dim1)]

  b <- spobj[["coefficients"]][[as.character(call$covariate)]]
  r <- summary(spobj)$rho
  p <- spobj[["coefficients"]][[as.character(call$tlag)]]


  direct_sr <- round(median(diag(solve(diag(dim1) - r*W)*b)), digits=5)
  direct_sr_l <- round(quantile(diag(solve(diag(dim1) - r*W)*b), probs = 0.025), digits=5)
  direct_sr_h <- round(quantile(diag(solve(diag(dim1) - r*W)*b), probs = 0.975), digits=5)
  
  indirect_sr <- round(median(colSums((solve(diag(dim1) - r*W)*b)) - diag(solve(diag(dim1) - r*W)*b)), digits=5)
  indirect_sr_l <- round(quantile(colSums((solve(diag(dim1) - r*W)*b)) - diag(solve(diag(dim1) - r*W)*b), probs = 0.025), digits=5)
  indirect_sr_h <- round(quantile(colSums((solve(diag(dim1) - r*W)*b)) - diag(solve(diag(dim1) - r*W)*b), probs = 0.975), digits=5)
  
  total_sr <- direct_sr + indirect_sr
  total_sr_l <- direct_sr_l + indirect_sr_l
  total_sr_h <- direct_sr_h + indirect_sr_h
  

  direct_lr <- round(median(diag(solve(diag(dim1) - r*W)*(b/(1-p)))), digits=5)
  direct_lr_l <- round(quantile(diag(solve(diag(dim1) - r*W)*(b/(1-p))), probs = 0.025), digits=5)
  direct_lr_h <- round(quantile(diag(solve(diag(dim1) - r*W)*(b/(1-p))), probs = 0.975), digits=5)
  
  indirect_lr <- round(median(colSums((solve(diag(dim1) - r*W)*(b/(1-p)))) - diag(solve(diag(dim1) - r*W)*(b/(1-p)))), digits=5)
  indirect_lr_l <- round(quantile(colSums((solve(diag(dim1) - r*W)*(b/(1-p)))) - diag(solve(diag(dim1) - r*W)*(b/(1-p))), probs = 0.025), digits=5)
  indirect_lr_h <- round(quantile(colSums((solve(diag(dim1) - r*W)*(b/(1-p)))) - diag(solve(diag(dim1) - r*W)*(b/(1-p))), probs = 0.025), digits=5)
  
  total_lr <- direct_lr + indirect_lr
  total_lr_l <- direct_lr_l + indirect_lr_l
  total_lr_h <- direct_lr_h + indirect_lr_h

matcell <- c("Direct Short-Run:", "Indirect Short-Run:", "Total Short-Run:", "Direct Long-Run:", "Indirect Long-Run:", "Total Long-Run:",
             direct_sr,  indirect_sr, total_sr, direct_lr, indirect_lr,total_lr,
             direct_sr_l, indirect_sr_l, total_sr_l,  direct_lr_l,  indirect_lr_l,  total_lr_l,
             direct_sr_h, indirect_sr_h, total_sr_h,  direct_lr_h,  indirect_lr_h,  total_lr_h)
fx <- as.data.frame(matrix(matcell,nrow=6,ncol=4))
colnames(fx)<-c("Effect", "Median Size", "2.5% Lower CI", "97.5% Higher CI")
  return(fx)

}




slx_plot_knb <- function(modResults, var){
       require(stringi)
       library(ggplot2)
       library(dplyr)
       
       modSumm=lapply(modResults, 
                      function(x) FUN=summary(x)$coefficients[,c('Estimate','Std. Error','Pr(>|t|)')])
       est <- array(NA, c(length(modSumm), 3))
       
       for (i in 1:length(modSumm)){
              
              est[i, 2:3] <- modSumm[[i]][var,c('Estimate','Std. Error')]
              est[i, 1] <- i +1
              
       }
       
       est <- as.data.frame(est)
       colnames(est) <- c("knb", "beta", "se")
       est <- est %>%
              dplyr::mutate(ad = .5/se) %>%
              dplyr::mutate(beta = beta*ad,
                            se = se*ad) %>%
              dplyr::mutate(lo95 = beta - 1.96*se,
                            hi95 = beta +1.96*se,
                            lo90 = beta - 1.645*se,
                            hi90 = beta +1.645*se) %>%
              dplyr::mutate(color_sig = ifelse(lo95 < 0 & hi95 > 0, "#ef8a62", "#2c7bb6"))
       
       p <- ggplot(est, aes(x=knb, y=beta, color=color_sig)) +
              geom_hline(aes(yintercept=0), linetype=2, color = "black") + 
              geom_point(size=4) + 
              geom_linerange(aes(ymin=lo90, ymax=hi90),alpha = 1, size = 1.5) + 
              geom_linerange(aes(ymin=lo95,ymax=hi95),alpha = 1, size = .5) + 
              theme_bw() +
              scale_x_continuous(breaks = seq(min(est$knb), max(est$knb), 1)) + 
              theme(legend.position="none",
                    legend.title=element_blank(),
                    axis.text = element_text(size=14),
                    text = element_text(size=14),
                    plot.title = element_text(hjust = .5, size = 16, face = "bold"))
       return(p)
       
}

